
!------ calculation parametrs of triangle (coef of basic functions) 

! THE area(n,it) is inverse of "real" DOUBLE area (=1./det) =>
! area (n,it) ( in program) = 1./(2 real_area)
! real_area = 0.5* (1./area(n,t))

subroutine init_areas
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'


do 13 i = 1,nx-1
    do 13 j = 1,nz-1

        ! Coordinates
        x1 = cord (j  ,i  ,1)
        y1 = cord (j  ,i  ,2)
        x2 = cord (j+1,i  ,1)
        y2 = cord (j+1,i  ,2)
        x3 = cord (j  ,i+1,1)
        y3 = cord (j  ,i+1,2)
        x4 = cord (j+1,i+1,1)
        y4 = cord (j+1,i+1,2)

        ! (1) Element A:
        det=((x2*y3-y2*x3)-(x1*y3-y1*x3)+(x1*y2-y1*x2))
        det1 = 1./det
        area(1,j,i) = det1

        ! (2) Element B:
        det=((x2*y4-y2*x4)-(x3*y4-y3*x4)+(x3*y2-y3*x2))
        det1 = 1./det
        area(2,j,i) = det1

        ! (3) Element C:
        det=((x2*y4-y2*x4)-(x1*y4-y1*x4)+(x1*y2-y1*x2))
        det1 = 1./det
        area(3,j,i) = det1

        ! (4) Element D:
        det=((x4*y3-y4*x3)-(x1*y3-y1*x3)+(x1*y4-y1*x4))
        det1 = 1./det
        area(4,j,i) = det1
        if( area(1,j,i).lt.0.or.area(2,j,i).lt.0.or.area(3,j,i).lt.0.or.area(4,j,i).lt.0 ) then
            call SysMsg('INIT_AREAS: Negative area!')
            write(*,*) 'x1,y1 = ',x1,y1
            write(*,*) 'x2,y2 = ',x2,y2
            write(*,*) 'x3,y3 = ',x3,y3
            write(*,*) 'x4,y4 = ',x4,y4
            write(*,*) 'area1 = ',area(1,j,i),j,i
            write(*,*) 'area2 = ',area(2,j,i),j,i
            write(*,*) 'area3 = ',area(3,j,i),j,i
            write(*,*) 'area4 = ',area(4,j,i),j,i

            stop 41
        endif
13 continue

return
end

!  1 - 3
!  |   |
!  2 - 4
!
!  diagonal / :
!
!   A:        B:
!
!  1---3         1
!  | /         / |
!  2         2---3
!
!  diagonal \ :
!
!   C:        D:
!
!  1          1---3
!  | \         \  |
!  2---3          2


function total_area( iph )
include 'precision.inc'
include 'params.inc'
include 'arrays.inc'

area_t = 0
do i = 1,nx-1
    do j = 1,nz-1

        if( iph .ne. 0 ) then
            if( iphase(i,j,phasez(j,i)) .ne. iph ) cycle
        endif
        x1 = cord (j  ,i  ,1)
        y1 = cord (j  ,i  ,2)
        x2 = cord (j+1,i  ,1)
        y2 = cord (j+1,i  ,2)
        x3 = cord (j  ,i+1,1)
        y3 = cord (j  ,i+1,2)
        x4 = cord (j+1,i+1,1)
        y4 = cord (j+1,i+1,2)

        ! (1) Element A:
        det=((x2*y3-y2*x3)-(x1*y3-y1*x3)+(x1*y2-y1*x2))
        area_t = area_t + det/2

        ! (2) Element B:
        det=((x2*y4-y2*x4)-(x3*y4-y3*x4)+(x3*y2-y3*x2))
        area_t = area_t + det/2

    end do
end do

total_area = area_t

return
end
